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DEVELOPMENTS IN THE COMPUTATION OF TURBULENT BOUNDARY LAYERS 
Morris W. Rubes In 

Ames Research Center, NASA, Moffett Field, California 94035, U.S.A. 


SUMMARY 

Computational techniques applicable to turbulent boundary layers are classified into solutions of 
Reynolds-averaged equations, in which all the effects of the turbulence are modelled, and solutions of three- 
dimensional, time-dependent Navier-Stokes equations, in which the large eddies are calculated and only the 
turbulence at scales smaller than the computational mesh spacings has to be modelled. Current computation 
costs place engineering computations in the first of these categories; large eddy simulations are aDpropriate 
currently for special studies of the dynamical processes of turbulence in idealized flow fields. It is shown 
that the two methods are interrelated and that each can gain from advances in the other. The degree of suc- 
cess of a pair of increasingly complex Reynolds stress models to broaden their range of applicability is 
examined through comparisons with experimental data for a variety of flow conditions. An example of a large- 
eddy simulation is presented, compared with experimental results, and used to evaluate the models for pres- 
sure rate-of-strain correlations and dissipation in the Reynolds-averaged equations. 


NOMENCLATURE 


a 

transverse body radius 

V 

subgrid fluctuating velocity component in ith 
direction 

Ci 

modelling coefficient in pressure rate-of- 
strain correlation, turbulence-turbulence 
interaction (Rotta term) 

u" 

fluctuating velocity in mass-weighted variables 



v w 

surface mass- transfer normal velocity 

C f 

skin- friction coefficient 




v' 

fluctuating velocity normal to surface 

D 1j 

turbulence production tensor, Eq. (37) 

y 

distance normal to surface 

c 

turbulence kinetic energy 

a 

scaling factor in turbulence simulation 

h 

static enthalpy 




a 

modelling coefficient in pressure rate of strain 

h M 

fluctuating entl alpy in mass-weighted 
variables 


correlation turbulence mean-fow interaction 



8 

modelling coefficient, Eq. (31), or scaling 

i 

length scale 


factor in turbulence simulation 

*m 

mixing length 

8* 

modelling coefficient, Eq. (31) 

P U 

turbulence production tensor, Eq. (37) 

6 

modelling coefficient in pressure rate of strain 
correlation - turbulence mean- flow interaction 

Pr L 

Prandtl number for molecular motions 



L 


Y 

modelling coefficient, Eq. (31) 

p* 

mean specific static pressure 

Y* 

modelling coefficient, Ea. (31) 

p 

static pressure 

Y 

modelling coefficient in pressure rate of strain 

p* 

specific static pressure, p/o 


correlation - turbulence mean-flow interaction 

p* 

resolvable fluctuating specific static 

6 

boundary-layer thickness 


pressure 

s u 

Kronecker delta 

p* 1 

subgrid fluctuating specific static 



pressure 

c 

eddy viscosity (turbulence kinematic viscosity) 

Q j 

heat flux vector 

c 

turbulence kinetic energy dissipation rate 

q 

turbulence speed 

E 1j 

dissipation rates of Reynolds stress component 

Re 

modelling coefficient, near-wall modifica- 


T U 


tion, Eq. (31) 

\ 

boundary- layer edge turbulence intensity factor 

Re T 

turbulence Reynolds number, Eq. (32) 

X 

low Reynolds number modelling coefficient. 

1 


Eq. (31) 

R u 

modelling coefficient, near wall modifica- 




tion, Eq. (31) 

u 

fluid viscosity 

S 1J 

rate of strain tensor 

V 

fluid kinematic viscosity 

U 1 

mean velocity component in ith direction 

v eff 

effective eddy viscosity in subgrid model 

U 1 

velocity component in the ith direction, 
mass-weighted in compressible flows 

p 

fluid density 


p 

instantaneous fluid density 

Q 1 

resolvable fluctuating velocity component 
in 1th direction 

0 

modelling coefficient, Eq. (31) 


REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 
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o* modelling coefficient. Eg. (31) 

Reynolds stress component 
vortlclty tensor 

u turbulence specific dissipation rate 

Subscripts: 

e boundary- layer edge condition or experlmen' 

tal value 

w surface quantity 

1,2,3 Cartesian axes direction 


( ) , a partial derivative with respect to the 
’ 1th coordinate 

Superscripts: 

T total quantity, turbulent plus molecular 

process 

(‘) partial derivative with respect to time 

Other: 

< > filtered average 

(*) time averaged 


INTRODUCTION 

Advances In computer technology and numerical analysis during the past decade have made it possible to 
compute the characteristics of turbulent flow fields with a degree of detail that was impossible in the past. 
This computational power has been applied to problems both In engineering and in basic fluid mechanics. 
Engineering methods have been confined largely to the solution of statistical equations of turbulence, 
usually for steady-state conditions. The Increased computational power has permitted the use of second-order 
closure methods wherein partial differential equations are used to describe the scales, intensity, and even 
the individual components of Reynolds stresses distributed throughout the flow field. It is becoming quite 
standard in advanced engineering problems to use two-equation models representing the transport of turbulence 
kinetic energy and a measure of the turbulence scale to establish the local eddy viscosity. The underlying 
Impetus to this work has been the premise that the increased complexity of a model tends to broaden its range 
of applicability, thereby making it a predictive tool. The past simple models of statistical turbulence, 
such as the mixing-length models, largely were used to explain flow-field behavior after the experimental 
results were obtained; they could be used only to interpolate between or moderately extrapolate conditions 
of a particular experiment. New situations required new experiments to guide modelling changes. With the 
more detailed second-order closure models, however, the rather large number of experimental coefficients 
employed requires drawing upoi experiments of different kinds of flow fields for evaluating the modelling 
coefficients of the different mechanisms. For example, the coefficient for the dissipation of turbulence 
kinetic energy is determined in part from experiments dealing with the decay of isotropic turbulence. Coeffi- 
cients for terms representing the exchange between Individual components of Reynolds stress by the correla- 
tion of fluctuations in pressure and the instantaneous rate of strain come largely from experiments in 
homogeneous turbulence created by either uniform shearing or normal strains. Thus, when these models are 
used in boundary-layer flows, many of their terms reflect the behavior of turbulence under other conditions, 
thereby possibly broadening the model's range of applicability. On the other hand, this reliance on a group 
of different kinds of experiments to establish the modelling coefficients often results in a somewhat less 
accurate representation of a particular flow field than is provided by a fine-tuned, simple empirical model. 
Engineers working continuously with certain kinds of flow fields tend to fine-tune the second-order models 
as well, without (it is hoped) losing too much of the generality potential within the model. 

Perhaps an even more important application of the powerful computation tools available today has been in 
the rather new field of the numerical simulation of the large eddies of turbulence. In these calculations, 
the three-dimensional and time-dependent character of the turbulent flow fields is retained. The principal 
approximations involved In these methods is In the manner of accounting for the scales of turbulence that are 
too small to be resolved, even in the largest of the computers, for the time-dependent and spatial ly-dependent 
boundary conditions, and for the initial field of the turbulence. Because the Initial and boundary conditions 
Involve so many degrees of freedom, it cannot be expected that individual computational realizations will be 
significant or even realistic. What can be expected or at least hoped for is that the results of the compu- 
tations viewed statistically will accurately reflect the highly nonlinear mechanisms that govern the dynamics 
of the flows. 

The computations, to date, have shown much promise but they need considerable development; moreover, they 
are much too costly to be considered as engineering tools. They are invaluable, however, as a technique for 
the study of fluid mechanics in that they yield a mass of information about a flow field that experiments 
involving discrete numbers of probes cannot possibly provide. The numerical analyst — when faced with making 
the choices necessary for starting his problem, the ranges of scales he is to examine, and the techniques of 
accounting f or the subgrid scales - is forced to consider the details of turbulence that modelers of statisti- 
cal turbulence have had to ignore. The apparent turbulence that numerical instabilities or bifurcations can 
create forces the serious worker in turbulence simulation to continually compare his numerical results with 
experimental data for similar flow fields. Because of the mass of detail contained within the calculations, 
where for example any statistical moment can be generated, It is often found that even the classic experiments 
lack data that unambiguously define the turbulence that was present. 

The purpose of this paper is to demonstrate that both methods of turbulence computation have much in 
coimion and that they are distinguished primarily by the fraction of the turbulence that is chosen to be 
modelled. The statistical methods model all of the turbulence. Ignoring most of the scale and all of the 
phase character of the actual turbulence. The large eddy simulations compute the actual physical character 
of the larger scales of the turbulence and model only those scales of turbulence smaller than the computa- 
tional mesh dimensions. The fraction of turbulence the subgrid model represents depends largely on the local 
turbulence Reynolds number and on the number of mesh points that the computer ran handle. 

In the present discussion, the two methods are Interrelated; a review is given of the success or failure 
of a pair of second-order closure methods to model the statistical properties of a variety of turbulent flow 
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fields, without adjustment of modelling constants; and some examples of large eddy simulations for simply 
strained homogeneous turbulent flow fields are presented and compared with data. Examples of statistical 
model Information that can be gained from these computations is shown. Finally, a plea is made for coordi- 
nated experiments and large eddy simulations, which together should prove to be most valuable in explaining 
the physics of turbulence under a variety of flow conditions. 

TURBULENT FLOW EQUATIONS 

In the analysis of turbulent motions. It Is generally believed that the basic physics of the fluid flow 

is contained within the Navier Stokes equations. Since these same equations apply at a point in space and 

time for both laminar and turbulent motions, the distinction between these types of flows arises from the 
Initial and boundary conditions the flows experience and from their response to small disturbances that are 
always present in real flow fields. This response is largely dependent on the Reynolds number of the flow. 

The properties of turbulence, then, are the consequence of the fluid instabilities that occur at high Reynolds 
number and the subsequent nonlinear, apparently chaotic, mixing processes that take place. It is these non- 
linear processes that produce a broad range of length scales of motion within the turbulent flow and this, in 
turn, affects finite difference computations greatly. 

The length scales range from those comparable to the characteristic dimensions of the apparatus down to 
those where the turbulent motions have largely been dissipated by viscosity into heat. Even the largest 
available computers fall far short of being able to resolve suer a broad range of scales for flow fields of 
technological interest. Although the prospects of increasing fie resolution of the turbulence scale with 
future computers is good (Ref. 1), it is not expected that it will be possible to compute the smallest dissi- 
pation scales, at Reynolds numbers of interest, in the reasonably near future. 

A variety of turbulence models has been developed to account for these small Irresolvable subgrid scr’es. 
These models have a great deal in common with the models for Reynolds stresses in statistical turbulence 
theory. This is demonstrated in this section through the equations for an incompressible fluid that describe 
the small scales of turbulence and their effect on the larger scales and the mean flow. 

The Instantaneous motion of an incompressible, viscous fluid Is described by the continuity and Navier- 
Stokes equations 


u j,J * 0 (1) 

and 

*1 * < u i“j ♦ «ijP* - ’ Mtt 1.j).j * 0 (2) 

where 

P* * p/o (3) 

The Instantaneous, local velocity can be expressed as the sum of three components: the time mean velocity; 
the sum of the fluctuating turbulence components whose length scales can be resolved by the finite-difference 
computational scheme; and the sum of those fluctuations too small to be resolved, namely 

U 1 “ U i + “i 4 uj (4) 


The other dependent variable, the pressure, can also be resolved in a similar manner 

p* = p* + p* + p*' (5) 


To convert Eqs. (1) and (2) to contain only resolvable dependent variables, it is necessary to average them 
or filter them in some manner. For the purposes at hand, it is not necessary to define the filtering process 
precisely. It can represent a weighted average over a line in space, a surface, a volume, or even a charac- 
teristic time comparable to the time scales of the small length srales of turbulence. The filtering process, 
represented by the symbol < >, will be defined in such a way as to accomplish the following: 

<uj> * 0 | 


<uJUj> = 0 
<ujUj> » 0 


( 6 ) 


<u 1.j > - <tt 1 > .j f 

Note that the fourth and fifth definitions in (6) imply that the two scales of turbulence are uncorrelated 
over the filter domain and that the domain is small compared to the scales of the mean motion. Leonard 
(Ref. 2) has demonstrated the limitations of these assumptions, but the simplicity resulting from them is 
attractive for the present development. 


j.j * U j.J T u j.J 


♦ ui 


0 


With Eq. (4), Eq. (1) becomes 


(7) 


( 8 ) 
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» 

When filtered according to Eq. (6), Eq. (7) reduces to 


U J.J + fi j.j * ° 

When Eq. (8) 1$ time-averaged In the Reynolds sense, that Is, 



f * 11m i 

f f(t)dt 

(9) 


T“*» 4 



there results 


u j.j 

“ t 

(10) 

Of course, subtracting Eq. 

(10) from Eq. (8), and Eq. 

(8) from Eq. (7) yields 



“J.J 

* 0 

(ID 

and 


U j.J 

■ 0 

(12) 


Similar operations of filtering and time-averaging on the Navler-Stokes equations (Eq. (2)) yield the fol- 
lowing forms of the momentum equation. The mean-momentum equations are 


Vu 


-Vtj + vU i.j,j - <vj 


(13) 


The momentum equations for the resolvable turbulence scales needed to evaluate the u^ in Eq. (13) are 

6 i + Vi , j " - + ' Vi.j ' ( Vj - ¥?.j - (<u i u j > - <_ W ) ,j (14) 


The equations for the Instantaneous values of the filtered moments of the subgrid scales, themselves at 
resolvable scales, can be expressed as 

<M i“k > + (u j + “j )<u i u k > ,j * - <u k u j >(u i + G i>.j - <u i u j >(u k + V,j 


* <u i u k u j > ,J * <u k p !i > * <u i p tk > 

* v<u i u k\j,j ‘ 2v<u i.J u i.J > (15) 

An equation for the subgrid turbulence kinetic energy, defined as 

<q 2 > « <u|u|> (16) 

follows from the trace of Eq. (15): 

<q 2 > + (Uj + Uj)<q 2 > j * -2<ujuj>(U 1 ♦ flj)j - <Uj(2p*’ + q 2 )>j 

* v<qa> .J.J “ 2v<u i.J u i.j > (17) 


Equation (13) shows that the influence on the mean motion of the two scales of turbulence Is through the 
sum of two Reynolds stresses, each associated with the different scales. In statistical turbulence theory, 
the different scales of turbulence are ignored by summing these Reynolds stresses into a single stress which 
is then modelled. Thus, all the effects of the turbulence on the mean motion are modelled. In large eddy 
simulations, however, the u^ are calculated as functions of time and in three dimensions and the correspond- 
ing large eddy Reynolds stresses are then computed through time averaging. (Actually, most large eddy simu- 
lations compute the sum of and u^, but in principle the mean flow is affected as stated.) Only the small 
scales are modelled and their influence is felt on the mean flow through the Reynolds stresses they contribute 
and in their effect on the larger scales of the turbulence. As the larger fraction of the turbulence spectrum 
Is computed, less reliance has to be placed on contributions of the turbulence model. A mechanism exists, 
then, for converging on the correct statistical description of turbulence through a systematic increase in the 
fraction of the resolved turbulence scales. It is not clear, to date, how the quality of the subgrid turbu- 
lence model affects this convergence rate and actually how many scales of the large eddies require computation 
to provide a good description of the turbulence transport mechanism in a variety of flow fields. Another 
apparent advantage of the method of large eddy simulations is based on the optimism regarding the generation 
of a universal subgrid model. This optimism reflects the experimental evidence that the small scales of tur- 
bulence In a variety of flow fields exhibit similar spectral characteristics when scaled in the Kolmogorov 
sense. 


Equation (14) illustrates the Influence of the mean flow and the subgrid turbulence scales on the resolv- 
able scales of the turbulenci . It shows that the growth of the resolvable turbulence along a mean streamline 
is acted on by the turbulence pressure and viscous diffusion in the same manner that the corresponding terms 
act in the rean flow. The additional last three terms represent the Interactions between the resolvable tur- 
bulence and the mean flow, the components of the resolvable turbulence, and the subgrid scales of the turbu- 
lence. It Is the latter terms that must be modelled to close the calculation of the resolvable scales. 
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Incidentally, these calculations must be performed in three dimensions and in time for each component of 
the resolvable turbulence; it is a rather costly process on today's computers. 

1 uations (15) end (17) provide insight into the modelling of the subgrid turbulence scales needed to 
"close" Eqs (13) and (14). When either Eq. (15) or Eq. (17) for the subgrid scale moments is compared with 
the equations representing the Reynolds stresses or kinetic energy in a statistical turbulence formulation. 

It is noted that they have essentially the same form except that the mean flow in the statistical equations 
has been replaced by the Instantaneous large-scale motions, and the subgrid-scale moments have a time 
dependence. This strong similarity between the subgrid-moment equations and the Reynolds stress equations 
suggests that much of the experience gained with statistical modelling procedures eventually will be able 
to be applied to subgrid modelling. At present, the limitations of computer storage encourage use of the 
simplest of subgrid models, analogous to the first-order closure methods, such as constant eddy diffusivities 
or mixing-length models in statistical Reynolds stress methods. Equations (15) and (16), however, suggest 
that second-order closure methods applied in statistical methods over the past decade will have a role in 
subgrid closure as well. Computer limitations and costs, also, will restrict large-eddy simulation in the 
near future to simple flow fields. The results, however, will provide considerable insight into the physics 
of turbulence and will contribute to modelling of the statistical equations. Some preliminary studies of the 
latter are given later in this paper. 

STATISTICAL TURBULENCE MODELLING 

Since the 1968 Stanford conference on the computation of turbulent boundary layers (Ref. 3), statistical 
turbulence modelling for engineering applications has gone in two directions. The first has been the fine 
tuning of first-order closure methods, involving algebraic mixing-length models. This was accomplished by 
fitting the models to well-defined experiments with attached boundary layers experiencing pressure gradients 
and/or surface mass transfer. References 4 and 5 are examples of this approach. Although these methods yield 
excellent representations of the data within their range of application, their abilities to extrapolate beyond 
the ranges of the experiments that form their basis is questionable. Any introduction of additional length 
scales into the boundary- layer characteristics, such as a transverse radius of curvature comparable to the 
boundary- layer thickness or an injection slot dimension, requires considerable remodelling of the length 
scales. In the absence of experimental guidance, this remodelling has to be based on ad hoc assumptions. 
Further, in flow fields where changes in the mean flow are rapid, the assumption inherent in most first-order 
ck-ure methods, that the turbulence remains In equilibrium with the mean flow, may not be true. The recogni- 
tion of these limitations of first-order closure methods has led to considerable work in the second direction, 
namely second-order closure methods, where one or more of the characteristics of turbulence is represented by 
•1 partial differential transport equation. These methods are based largely on the concepts presented in the 
pioneering papers by Kolmogorov (Ref. 6), Chou (Ref. 7), and Rotta (Refs. 8,9). As explained in the Intro- 
duction, the impetus behind the development of second-order closure methods was the belief that they have the 
potential of a broad range of applicability and may, with further development, become predicti ve tools in 
engineering computations. 

An ideal predictive turbulence model would be one that could remain unaltered in form and in its empiri- 
cal coefficients for all flow fields. It is questionable, however, that such an ideal universal model can be 
achieved within the framework of statistical models, even when allowance is made for acceptable engineering 
error. Such models inherently ignore spectral and phase relationships between eddies of different sizes. 

The larger eddies in a turbulent flow are known to reflect the particular nature of the flow and this alone 
is sufficient to raise doubts regarding the potential universality of statistical models. Of course, zonal 
turbulence models that differ from flow to flow but can be related to some particular mean flow feature are 
also of value to the design engineer; they may be the best that can be exoected of statistical models. Such 
zonal models, however, introduce mathematical difficulties in the identification of the bounds of different 
zones of applicability of the individual models, and the means of coupling the interaction between these 
zones. It is much easier to deal with the same model throughout the flow. In view of this, it would be 
illuminating to learn how well or how poorly a fl xed model could work on a variety of boundary-layer flows. 

To demonstrate this, calculations based on a pair of fixed second-order closure models will be compared here 
with data from a variety of experiments. 

The particular models chosen for this comparison have been developed for the most part by Wilcox and 
Traci (Ref. 10) with some collaboration by the present author (Ref. 11), and were an outgrowth of the early 
work of Saffman (Ref. 12). One model uses an eddy viscosity, which is dependent on the kinetic energy of 
turbulence and the dissipation rate per unit of kinetic energy (a specific dissipation rate or Saffman's 
"pseudovorticity"). The other model closes the Reynolds stress equations directly, with the scale of turbu- 
lence again being defined with the specific dissipation rate. In choosing these models for comoarison here, 
the author does not wish to imply that he believes them to be the best of the second-order closure methods 
available today to represent boundary- layer flows. In some respects, the models developed by Launder and 
his colleagues (Refs. 13,14) are more general and represent the mean flow very close to a surface in a more 
realistic manner (Refs. 14,15). On the other hand, computations of some compressible flow fields with two- 
equation models favor the Wilcox-Rubesin model (Ref. 11). It is not clear at present which, if any, existing 
model will be most uniformly valid for all applications. The primary reason for presenting models in which 
the author was involved was his access to computer codes containing them; as a result, the cod’s could be 
used to generate the examples that follow. Also, since both of these models and those identified with 
Launder utilize essentially the same data to establish their modelling coefficients, it should not make much 
difference in the examination of the universality of second-order closure modelling of boundary-layer flows 
which family of models is employed. 

Mean Flow Equations and Boundary Conditions 

The mean flow equations used in computing the examples that will be treated later are written for a 
compressible fluid in mass-weighted-average dependent variables (Ref. 16). The conservation equations for 
mass, momentum, and energy are as follows: 


p + (oUj ) , j * 0 


(18) 
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(p*u 1 ) ♦ (pUjU^.j ■ -p, 1 ♦ ( pT 1j)*j ( 19 ) 

(ph) ♦ (pUjh),j * p + UjP.j + pxj^ - (pQj T ),j (20) 

Here, the symbols tIj and Qj T denote the specific mass-weighted-average total shear stress and heat flux 
that Include the contributions of both the molecular and turbulent transport. These quantities are defined as 


T 1J * 2v S 1j * 5 u k,k 4 1J * T 1J (21) 

and 

°j T * ' h *j + Q j (22) 

the mean rate of strain tensor In Eq. (21) Is 

S 1J * \ (u 1.J + “j.1* (23) 

Finally, and Qj are the mass-welghted-averaged Reynolds stress tensor and heat flux vector defined by 

<PU|"u "> 

T u - - — r 1 - 

<pu 4 "h"> 

where p Is the Instantaneous density, <e> denotes the time average of e, and e" Is the fluctuating part 
of e In mass-weighted-average formulation. The surface boundary conditions for Eqs. (18) to (20) are: 



at x 2 ■ 0 

uj * 0 

u 2 * 0 or v w (x t ) 

h * h w (xi) or ( 3h/ 3 x 2 ) * (3h/3x 2 ) w 


(25) 


All flow variables approach free-stream flow conditions In general flow-field computations. For the special 
case of two-dimensional boundary layers, boundary conditions at the boundary-layer edge are 


at x 2 « «(xi) 


uj * U e (x)) 


h • h e (x] ) 


(26) 


The two models used to close these equations are given in the following sections. 

Two-Equation Eddy-DI ffusi vity Model 

The two-equation model considered here utilizes an eddy diffusivity defined as 

c * y*e/u (27) 

where the turbulence kinetic energy e and specific dissipation rate u are given y the turbulence model 
ling equations: 


(p‘e) ♦ (pUje),j - p-r^jUj j - a*pu>e ♦ [ ( w + o*pc)ekj],j 


(28) 


and 


(p’u) 2 ) + (pUjU 2 ) , j * y p t 1 jU 1 j - [a + 2o(t, k t,| ( )]pu 3 ♦ C(u + Opt) (u 2 ) ,j] , j 


(29) 


To account for compresslbllltv, all the dependent variables are expressed as mass-weighted averages (Ref. 16). 
The length scale Is represented by 

e>/ 2 


t • 


(30) 
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The modelling closure coefficients employed are as follows 

B * 3/20, B* * 9/100, o*o*= 1/2 
y* * [1 - (1 - x 2 )exp(-Re T /R e )] 

yy* * y„[1 - 0 - x 2 )exp(-Re T /RJ] 
y. * 10/9, X - 1/11, R e - 1, R u - 2 
The Reynolds number of turbulence Is given by 


(31) 


Re 


T * 


e‘/ 2 t 

V 


(32) 


The boundary conditions appropriate to these modelling equations, when they are applied to boundary layers, 
have been guided by asymptotic analysis and reference to other models. The surface boundary conditions for 
Eqs. (19) and (20) are as follows: 


at X 2 * 0 


at x 2 = «(x,) 


e * 0 

20 v w 

10 -+ 

bx 2 2 ; 

e * iUg 2 (xj ) 
l * 0.09 6* 1 / 4 «(x 1 ) 


(33a) 


(33b) 


As the quantity 1 / B*'/ 4 behaves much like the classical mixing length, the proportionality coefficient of 
0.09 in Eq. (19) is read J ly seen to be consistent with the Escudier eddy-viscosity model (Ref. 17). 


Since It was desired to model Reynolds stresses that do not necessarily align with the mean rates of 
strain, the constitutive relationship relating these quantities was written as 

t u * I e6 u + 2e ( s u * i u k,k 6 u) + 1 7^r^s~s~) (Sin,n,nj + Sj,A,i) (34) 

“ + “rravW 

where the third term on the right was guided by the work of Saffman (Ref. 18). The vorticity tensor used 
here is defined as 

n u “i (u i,j * U j.i ) (35) 

and the mean rate of strain is given by Eq. (23). 


Reynolds Stress Equation Model 

The modelling in the Reynolds stress equations (RSE) presented here utilizes (1) a particular version 
of the pressure rate-of-strain correlation presented by Launder et al. in Ref. 14, (2) gradient diffusion 
for third-order correlations involving velocity and pressure, and (3) isotropic dissipation, following 
Launder et al., the pressure rate-of-strain correlation is represented as 


pTu 


l.j 


“j.l* 


C^uilT 


1j 


♦ ! 4 u e ) - 4 ( p 


1j 


! P5 u) 


K°ij ’ f P °i.i) * ;eS 


ij 


(36) 


where the first term on the right, called the Rotta term, is proportional to the anisotropy of the turbulence. 
The terms preceded by the modelling coefficients a, §, and y are contributed by the interaction of turbu- 
lence and mean flow expressed in terms of 


P 1j ” T 1k u j.k + T jk u i,k 
D 1j ' T 1k u k,j + T jk u k,1 

P * l P ij * I D 11 ■ T mn S nm 


(37) 


In choosing the values of a, 6, y» Wilcox and the author opted to rely on experimental data rather than the 
symmetry arguments recommended by Rotta (Ref. 8) and carried out by Launder et al . (Ref. 14). The first 
experimental observation employed was that 


r 33 s } < T n + T z 


(38) 
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when a homogeneous turbulent flow is equally stretched in the x 3 direction and compressed in the x 2 
direction, which by continuity in an incompressible fluid leaves the x 3 direction unstrained (Ref. 19). 

The second observation was that a field of homogeneous turbulence in rigid body rotation decays without 
developing anisotropy (Ref. 20). The latter data forces a = §. The remaining constants are evaluated, 
after representative values of tjj/t,, and ri 2 /e are established from data in a homogeneous shear flow or 
in the law-of-the-wal 1 region of a flat-plate boundary layer. The model used here was based on the approxi- 
mate relationships t 22 /tu 5 1/2 and r 12 /e * 0.3, both of which are consistent with the form of the two- 
equation model in shear flow. Further remarks regarding the modelling constants found in this way will be 
made in the section on large eddy simulations. 

With the modelling described, the Reynolds stress equations expressed in terms of the components of 
Reynolds stress are 


(ct 'lj ) * (p Vlj ) 'k * - pT 1m u j,m * pT jm u 1,m + ! 6V,e6 1j * C i B * p “( T 1j + f efi ij) 

+ p ( T jm S mi + T 1m S mJ * t T mn S nm 4 1j) + I p6 ( S 1j ' I u k,k 6 1j) + C(p + °* pe)T 1j,k ] -k 


(39) 


The components of the Reynolds heat flux are modelled with 

(°Q 1 ) + (pu j Q 1 ) »j * P T lj h *j - pQjUi.j - B**puQi + j(p^-+ o**pe)q l 

The specific dissipation rate used to provide the scale of the turbulence is again given by Eq. (20). 
In these equations 


’j 


(40) 


e * e/w 
e “ - ? T 11 


(41) 


and the modelling coefficients are 


where 


6 * TO • 6 * * TM * 6 ** “ 75 • 0 ’ °* * 7 • °** " 2 


*j * C u 1 - (1 - *2)exp^lj 


(1 - X 2 )exp^ 
(1 - x 2 )exp( 


CT ' t. f 1 * " ' 

c i. = [f - f exp(-5x)J 
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Again, the Reynolds number of turbulence Rej Is given by Eqs. (30) and (32). At the solid surface, 
Eqs. (24) again apply and, in addition 


V° ) 

T 1j - 0 I 

At the boundary- layer edge, in addition to Eqs. (24), it. is required that 


(43) 


at x 2 * 6 


v° 


T 1j * I ' u e 2 (*i)«ij 


(44) 


Examples of Turbulence Model Application 

As an initial example of the use of the turbulence models presented here, consideration is given to the 
distortion of a field of fully developed, homogeneous turbulence by application of plane or normal strains. 
This case acts as a test of the models when near-surface effects are absent. The particular case treated 
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corresponds to the flow in the experiment by Tucker and Reynolds (Ref. 19). Figure 1 shows a schematic 
diagram of the test channel. The fluid enters at the left, is conditioned through screens, and passes 
toward the right In parallel flow until it reaches the station where the constant rate of strain is applied. 
The constant strain is achieved by exponentially expanding the channel in the x direction and contracting 
it in the y direction so as to maintain a constant cross-sectional area. The straining causes the initially 
nearly isotropic turbulence to become anisotropic, a measure of which is the straining parameter plotted as 
the ordinate in the lower portion of the figure. After the fluid is strained, it is returned to a parallel 
flow. Here the fluid tends to return to isotropy. 

Use of the RSE model in computing this flow is noted to yield reasonably good agreement with the data. 

At the initiation of straining, the computed growth of anisotropy is somewhat faster than the data, but this 
trend reverses toward the end of the straining region and downstream. These trends were influenced by the 
assumed form of in Eq. (42) which was adjusted to fit an aggregate of homogeneous flow experiments, not 

just that of Tucker°°and Reynolds. Use of the two-equation model in the computations suffers in two respects. 
First, the model shows the difficulties of all eddy viscosity models when a sudden application or removal of 
mean strain occurs. Although the elements making up the eddy viscosity (e and u) vary continuously where 
the discontinuities in strain occur, the corresponding Reynolds stresses are still discontinuous. Second, 
the rate of generation of anisotropy, by the two-equation model away from the stations where the discontinui- 
ties occurred, was too slow. Thus, although the RSE model has been shown to yield rather good agreement with 
the data, this example illustrates that two-equation models are limited to flow conditions with more gradual 
application of mean strain. 

An example of the application of the turbulence models to a boundary- layer flow of an incompressible 
fluid is demonstrated in Figs. 2 through 4. The data are from the experiment conducted by Bradshaw in which 
a turbulent boundary layer was exposed to a sudden application of an adverse pressure gradient (Ref. 21). The 
data points designated with open symbols result from a reinterpretation of the basic data by an independent 
analysis (Ref. 22) and provide ar. indication of the uncertainty inherent in the data. Figure 2 shows the 
distribution of skin friction and boundary- layer shape factor along the test zone. The computations were 
started by matching the momentum-thickness Reynolds number from a flat-plate calculation to the Reynolds num- 
ber measured at the upstream station. Beyond this station, the experimental pressure distributions were 
imposed on the boundary- layer calculations. It is observed that the two-equation model and the Reynolds- 
stress model both yield skin friction and shape factor results that are nearly the same, and that both agree 
well with the data. It should be noted that no adjustments were made to the modelling to account for the 
pressure gradient. 


The measured and computed mean-velocity profiles at the farthest downstream station are shown in Fig. 3, 
which is expressed in wall-law coordinates. The computations based on both models yield results in good 
agreement with the data. In the "law-of-the-wal 1 " region, the computed results agree with the standard 
logarithmic formula. In that region, also, there is a little better agreement with Coles' reinterpretation 
of the data (Ref. 22). Neither of these observations is surprising, as the use of the logarithmic law with 
the constants shown played a major role in the data reinterpretation and in establishing some of the modelling 
coefficients used in both models. Perhaps more significant is that the computed results based on both models 
show the enhanced contribution of the "wake" region that is characteristic of boundary- layer flows in adverse 
pressure gradients. 

Figure 4 compares the Reynolds stress, turbulence energy, and mean-velocity profiles, computed from the 
two models, with Bradshaw's data. In these figures, the distance normal to the surface has been normalized 
by the boundary- layer thickness. Generally, the Reynolds stress components and the turbulent kinetic energy 
given by the two models differ less from each other than from the data. The normal components of Reynolds 
stress — <u' 2 >, <w' 2 > and, to a lesser degree, the kinetic energy — are evaluated rather ooorly in the inner 
half of the boundary layer. On the other hand, the normal Reynolds stress <v' 2 > and the shear stress 
<-u'v'> fit the data much better. This reflects the adjustment of some of the .odelling coefficients to 
provide good mean-velocity profiles in a flat-plate boundary layer, where a good evaluation of the Reynolds 
shear stress is paramount. In the coordinates of this figure, the two-equation model shows a little advan- 
tage over the RSE model; however, both models yield results that reflect the "flattening" of the velocity 
profile introduced by the adverse pressure gradient. 


The results that have been shown here are representative of comparisons of many other sets of data with 
computed results based on the two models. For attached, subsonic boundary layers on flat plates, with or 
without pressure gradients, there seems to be no advantage to the Reynolds stress model. Under these flow 
conditions, the departure of the turbulence from being in equilibrium with the mean flow is apparently too 
small to cause the two-equation model difficulties that were indicated earlier with suddenly distorted 
homogeneous flows. 


The growth of a turbulent boundary layer on a small-radius circular cylinder with its axis parallel to 
the free stream is an example of a flow where the mixing-length formulation required considerable change, from 
that on a flat plate, to make computation* conform with the experimental results. Generally, the data showed 
that as the ratio of the boundary- layer thickness to the transverse radius increased, the "wake" region of the 
boundary layer diminished (as occurs on a flat plate in a favorable pressure gradient). Also the skin-friction 
coefficient at a given Reynolds number based on boundary- layer thickness increased with diminishing body 
radius. When Rao (Ref. 23) examined this flow field, he concluded that conformance with his data could be 
achieved if he employed a wall-region law equivalent to setting the mixing length in the inner region of the 
boundary layer to 


[tn(i + n 

^/a)/l + y/a 

L < 

y/at J 


(45) 


In studying the same problem, Richmond (Ref. 24) deduced a formulation equivalent to 


i 


m 


ky 


0 ♦ y/2a) ' 

(1 + y/a) 3 / 2 
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For values of y * a, Eq. (45) remains within a few percent of t m = ky, whereas Eq. (46) shows a reduc- 
tion of about 50%. These differences result in adjusting the position of data in opposite directions, in 
terms of y + , relative to the universal "law of the wall" and demonstrates vividly how arbitrary the exten- 
sions of mixing-length closure can be. 

The RSE second-order closure model given here was applied without change to this type of flow to see if 
the experimental data could be represented. Figure 5 compares the computed results and measured data in terms 
of the effect on skin friction of the ratio of «/a. Four sets of data are utilized (Refs. 24-27) to cover 
a sizable «/a range. The computational results are shown as a band because the ordinate employed does not 
collapse all of the Reynolds number dependence. The computations with the unchanged RSE model represent the 
trends of the data well up to a value of «/a * 10, where the C* is 30% higher than that on a flat plate. 
Beyond this, a significant departure occurs from the data of Willmarth et al . (Ref. 27) where the measured 
Cf is increased to more than twice that of a flat plate. Apparently, the RSE model in its present form fails 
to fully account for changes in the wake character of the boundary layer over a body with an extremely small 
transverse radius. 

The next examp 1 a to be treated deals with the topic of streamwise curvature, the importance of which was 
first recognized by Bradshaw (Ref. 28). The Reynolds-stress equations were applied to this problem direct’y 
through the conversion of the coordinates from Cartesian to curvalinear, with one axis tangent to the surface, 
s, and the other normal to the surface, n. The two-equation model, however, required reinterpretation of the 
meaning of the symbol e, treated as the kinetic energy earlier. Details of these transformations are given 
in Ref. 11. 

For flow over a streamwise curved surface, the s, or curvilinear coordinate system, introduces terms in 
the Reynolds stress equations analogous to centrifugal and Coriolis forces in the momentum equations. When 
the normal stresses are added together, however, most of these additional terms cancel, resulting in an energy 
equation that is essentially the same as on a flat surface; the only change is in the production term where the 
mean-velocity gradient 3u/3y is replaced by (3u/3r) - (u/R). The specific energy-dissipation-rate equation 
is also changed in the same way. Thus, a direct application of the two-equation model as given earlier would 
not show streamline curvature effects of the magnitude indicated by a Reynolds stress model or by the experi- 
mental data (e.g.. Ref. 29). This deficiency was corrected by (1) observing that the Reynolds shear stress 
and v' 2 equations in the RSE model in s, n coordinates added similar terms because of the streamwise 
curvature and (2) identifying e with a "mixing energy" rather than a "kinetic energy." The symbol e then 
is redefined as 


e 



9 

T T nn 


(47) 


which follows from the basic model in a homogeneous shear flow where the turb ulenc e production and dissipation 
are in balance. With Eq. (47) and the Reynolds stress equation for v' 2 and u'v' as guides, the e equation 
for use with the two-equation model only is written in an ad hoc manner as 

ue, $ + ve, n + f ^ T * T ( u, n ‘ ft) * e * e “ + [ y + °*c)e» n ]. n (48) 

with 

t - c(u. n - £) (49) 

and, where c follows from Eq. (27); all the modelling coefficients and relationships employed in the two- 
equation model introduced earlier are retained. The third term on the left side of Eq. (48) represents the 
principal extra rate of turbulence production introduced by the longitudinal surface curvature. 

Calculations based on these model modifications for streamline surface curvature are compared in Fig. 6 
with data obtained by So and Mellor (Ref. 29) for a boundary layer on a convex wall with an adverse pressure 
gradient. The data represent the surface skin-friction coefficient and shape factor along the surface. The 
computations include the RSE and two-equation models, both with and without the corrections for longitudinal 
surface curvature. The computed results were matched to the first station by assuming the flow upstream of 
the station to be on a flat plate of a length to yield the co. rect skin friction there. The calculations 
with the RSE or two-equation model unmodified for streamwise curvature show little of the drop in skin- 
friction coefficient experienced in the experiment. The modified models, on the other hand, give an excellent 
representation of the skin-friction behavior. This is rather remarkable for the two-equation model, when 
its ad hoc formulation is considered. Finally, both modified models represent the form factor data also 
quite well . 

It may seem to be Illogical in the test of the universality of a turbulence model to make the modifica- 
tions indicated for introducing the effects of streanwise curvature. For the RSE model, the modifications 
were purely geometric and were introduced by selecting the appropriate coordinates for the problem considered. 
No physical modelling changes were made. The original two-equation model, on the other hand, was insensitive 
to changes in the coordinate system and an extra production term had to be added to the e equation. One 
can view the need for the change as an Indication of the inherent weakness of the original two-equation model, 
or the view of the final e equation as the basic boundary-layer model that then is simplified geometrically 
for planar surfaces. 

The remainder of the two-dimensional boundary layers considered here involve compressible flows. As 
noted earlier, the extension to compressible flows is achieved by adopting dependent variables that are Favre 
mass-weighted averages. In these variables, the conservation equations take forms that avoid terms that are 
explicitly dependent on the turbulent density fluctuations, Term-by-term, the equations are comparable to 
their Incompressible counterparts, with the compressibility entering only through the mean density variations. 
Although Favre averaged-model equations are, in the main, parallel term-by-term to their Incompressible 
counterparts, they have additional terms explicitly dependent on density fluctuations that require additional 
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modelling (Refs. 30,31). In practice, however, it was found that these additional terms could be ignored, 
even under conditions involving large pressure gradients where the terms are their largest (Ref. 32). 
Reflecting this, the models presented here neglect these additional terms. 

These models were applied to the calculation of skin friction on a cooled flat plate in airflow at a 
Mach number of 5. The results of these calculations based on the two models are compared with values given 
by the van Oriest II formulas in Fig. 7. Thes^ srmulas have been shown (Ref. 33) to represent the bulk of 
existing data under these conditions to about +10%. The agreement between all the methods is excellent, 
but this is not surprising in view of the similarity of the density scaling in Favre averaging (Ref. 11) and 
that which is inherent in the van Driest formulation. 

Compressible turbulent boundary layers experiencing severe pressure gradients are cases where the models 
are tested more severely. The first example of such a flow is the experiment conducted by Lewis et al. 

(Ref. 34) at a M?ch number of four. In that experiment, an axisymmetric jrbulent boundary layer on the 
adiabatic interior wall of a circular cylinder was subjected to an adverse pressure gradient followed by a 
favorable pressure gradient. The pressure gradients were achieved by means of a shaped center-body; a pres- 
sure rise of 9 times the upstream pressure was attained before pressure relaxation occurred. Figure 8 shows 
the distribution of the surface skin-friction coefficient within the test zone. The coefficient shown is 
defined in terms of the upstream boundary- layer edge conditions, not the local, and is therefore proportional 
to the surface shear. The Reynolds number at the initial station was about 7 * lo 6 . Along with the computed 
results based on the models presented here, computations based on the Marvin Shaeffer code (Ref. 35), which 
las been extended to contain a classic mixing-length model essentially identical to that of Cebe:i (Ref. 4), 
are given for comparison. The mixing-length model fails to capture the full rise of the skin friction 
caused by the adverse pressure gradient. On the other hand, it follows the data in the region of favorable 
pressure gradient quite well. The second-order closure models demonstrate a much better prediction of the 
rise in skin friction in the adverse pressure gradient region; however, in the following favorable pressure 
gradient region they show somewhat too large a drop in the skin friction. The two second-order models yield 
essentially equivalent results. 

Figure 9 shows data from a similar experiment conducted by Horstman et al. (Ref. 36) at an initial Mach 
number of 2.3 and over a large range of Reynolds numbers. In addition to computations based on the two models 
considered here, computed results from two versions of a mixing-length model and another two-equation model 
are also shown here. At this Mach number, M,„ = 2.3, the onset of an adverse pressure gradient first reduces 
the skin friction before a rise similar to that which occurred in Fig. 8 also occurs. Generally, the com- 
puted values of skin friction from all the models conform to the trends in the data caused by the change in 
Reynolds number and the effective pressure gradient. One exception is the behavior of the mixing-length 
model unmodified for pressure gradient which indicates separation at the lowest Reynolds number. At the higher 
Reynolds numbers, the difference between the modified and unmodified mixing-length model become very small. 

From this figure, conclusions regarding the relative merits of the different models would be indecisive. In 
Fig. 10, when the adverse pressure gradient is applied over a greater distance, the models behave in a somewhat 
different manner. At the lowest Reynolds number, the unmodified mixing-length model no longer indicates 
separation. Also, at the lower Reynolds numbers, the second-order closure models are in much better agreement 
with the data. Omission of the explicit density fluctuation terms resulting in the model equations after 
Favre averaging is justified by these examples. Incidentally, the computed results labeled Aeronautical 
Research Associates of Princeton (ARAP) are based on a Reynolds stress model utilizing primitive dependent 
variables Including the whrle gamut of fluctuating density terms (Ref. 37). 

The remaining examples of two-dimensional boundary layer and near-wake flows were computed with the 
compressible Navier-Stokes equations to account for strong interactions between the shear layers and the 
inviscid flow. Because these codes are costly to operate they have been limited, at least to date, to con- 
tain models of turbulence of the two-equation kind or simpler. Therefore, the RSE model will not appear in 
these examples. 

Figure 11 shows calculations of the surface pressure and skin-friction distributions compared with data 
in the region of the interaction of a normal shock wave with a fully established turbulent boundary layer 
(Ref. 38). The schematic diagram shows that the flow field was developed on the surface of a tube within a 
sligntly supersonic main flow. A normal shock wave was generated and positioned along the test section with 
a variable blockage device at the downstream end of the test section. The figure at the left shows that the 
computational results based on the two-equation model generally agree well with the measured surface pressure 
distributions at the five Mach numbers tested. The departures that exist from the data are small and incon- 
sistent enough to hide any systematic deficiencies in the computational model. The computed skin-friction 
coefficients again conform to the main features of the data. The calculations show a downstream movement of 
the minimum in skin- friction coefficient with increasing Mach number. If the extreme Mach number cases are 
emphasized, a similar movement is seen in the data, although of larger extent. The inaccuracies inherent in 
skin- friction measurements can possibly exaggerate the movement of the minimum skin-friction coefficient and 
could be the source of these differences. 

An example of a strong interaction between a boundary layer and a shock wave at higher Mach number is the 
experiment of Settles et al . (Ref. 39) with a turbulent boundary layer traversing a compression corner. The 
computations used for comparison with the data are from Ref. 40. Figure 12 shows computations and measure- 
ments of surface pressure and surface skin friction for two deflection angles of the compression corner, 
a = 20° and 24°. Besides the two-equation model under consideration, three other models have been used in 
these computations. Models not shown in any of the earlier examples are a kinetic energy model with an 
algebraic length scale (Ref. 31) and the Jones-Launder two-equation model (Ref. 13). A general observation 
is that both two-equation models yield essentially the same results, except for the level of skin friction in 
the reversed region. This suggests the kinds of measurement needed to distinguish between models. Comparison 
of the experimental data with the computations reveals that the two-equation models permitted the location of 
the onset of the increased surface pressure ahead of the compress 'on corner to be computed quite well. For 
a ■ 20°, these models yield excellent pressure distributions over _ne separated zone and on the deflected 
surface beyond reattachment. For a = 24°, the calculated pressure in the separated zone is somewhat high, 
although ahead of reparation and after reattachment the pressure is again evaluated quite well. 
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The zero-equation and one-equation models show late onsets of pressure rise, then overshoot the data, 
and then blend with the data far downstream. The comparison with the skin-friction data do not show all the 
same trends In the computations with the different models. The two-equation models define the position of 
the onset of the fall-off of skin friction, and the other models again lag this. The fall-off of the skin 
friction given by the two-equation models is faster than the data show, so that the points of separation are 
predicted upstream of where they actually are. The two-equation models also yield too long a separated 
region so that reattachment is calculated to be downstream of the experimental results. The comoarison is 
inconclusive, regarding which of the models best fits the downstream data, except that the algebraic model 
yields results that consistently fall lower than the data. 

Figure 13 compares the computations and measurements of the effect of Reynolds number on the extent of 
the upstream pressure influence ahead of the compression corner. The distances considered are shown schemat- 
ically in the left-hand sketch in Fig. 13. Two deflection angles are considered: a = 16° and a = 20°. 

It is seen that the observation made earlier that the two-equation models evaluate the position of the onset 
of the pressure rise best is borne out in the figure over the entire range of Reynolds number covered in 
the experiment. 

Another example in which >.omputat1ons with the two-equation model have been compared with experimental 
data is the work by Viswanath et al. (Ref. 40). The experiment was conducted at the trailing edge of a 
flat-plate test model that terminated with a 12.5° total included angle wedge. In the example cited here, 
the model and wedge trailing edge were both kept at zero angle of attack to an airstream at H = 0.7. 

Figure 14 compares measurements of the mean-velocity profiles just upstream and downstream of the trailing 
edge with computations employing the two-equation model and an algebraic model (Ref. 4) in both Navier-Stokes 
and boundary-layer equation. The high chord Reynolds number of 40 * 10 6 ensured a fully turbulent boundary 
layer well ahead of the trailing edge. In this figure, e 0 represents the momentum thickness of the boundary 
layer 0.4 cm upstream of the trailing edge; it is equal to 0.2 cm. The computations employing either model 
in the Navier-Stokes equations agree better with the data than the same models in the boundary- layer computa- 
tions. Under these conditions, either on the wedge or just beyond the trailing edge, the flow is more sensi- 
tive to the interaction between the shear flow and the inviscid flow regions than to the particular turbulence 
model. Apparently the rate of change of the mean motion, even this close to the trailing edge, is sufficiently 
slow for either an equilibrium-model or a two-equation model to still apply. Farther downstream in the wake, 
both models and both computation techniques yield essentially the same results. It is important to note that 
the sudden removal of a surface downstream of the trailing edge did not cause any difficulties with the near- 
wall modifications represented by Eqs. (31), as they blended smoothly toward their asymptotic values farther 
in the wake. 

The final example cited here is the response of a turbulent boundary layer to a sudden application of 
transverse shear, as studied experimentally on an axisymmetric rotating body in Refs. 25 and 26. A sketch of 
the model configuration is given in Fig. 15. The free-stream velocities in these experiments ranged from 
10 to 19 m/sec. A comparison of the data from the two experiments with a mixing-length model (Ref. 42) modi- 
fied with Eq. (45), the two-equation and RSE models, and the ARAP model (Ref. 37) has been represented in 
Ref. 43. In the computations with either of the eddy-viscosity models, it was necessary to introduce an 
additional assumption regarding the ratio of the eddy diffusivity corresponding to the transverse flow to 
that of the longitudinal flow. The need for assuming some value for the ratio is an inherent problem in the 
application of any scalar eddy-viscosity model to a three-dimensional boundary layer. In the computations 
with the two-equation model, this ratio was set equal to unity, as it was for the mixing-length model in 
Ref. 42. It was found in Ref. 43 that computations based on the simple mixing-length model yielded results 
in general agreement with the measurements of the mean flow. The two-equation and RSE models showed compari- 
sons that were only somewhat better than the simpler model. The improvement achieved by the second-order 
closure models seemed to be limited by too rapid a response to the transverse shear. The relative agreement 
between the RSE model and the scalar eddy- viscosity models can be explained by reference to Fig. 15 where the 
ratio of the eddy viscosities calculated from the RSE model are compared to the data from the two experiments. 
First, it is observed that the two similar experiments result in data in serious disagreement. The appro- 
priate ratio cannot be established experimentally. The RSE model, with or without the effects introduced by 
transverse curvature, shows that the ratio of eddy viscosities remains within ±105! of unity over most of the 
transverse boundary layer, as assumed in the eddy- viscosity models. Near the outermost edge of the transverse 
boundary layer, in the vicinity of the onset of the transverse shear, the ratio drops to a smaller number. 

As there is little momentum change near the boundary- layer edge, differences in eddy viscosity such as these 
have negligible effect on the transverse-velocity profiles. Thus, the choice of the eddy viscosity ratio of 
unity in the simpler models is not inconsistent with the evaluation of the RSE model. 

Concluding Remarks Regarding Statistical Modelling 

From the foregoing set of comparisons of e.verimental data and computations employing a pair of fixed 
second-order closure models, and with other models as well, it is observed that the second-order closure 
models generally have a broader range of application than do the algebraic closure models. The two-equation 
eddy- viscosity model is accurate over a large range of Reynolds numbers for attached two-dimensional incom- 
pressible or compressible boundary layers on impervious surfaces, even those with small zones of separation. 

The Reynolds stress model, in addition, has advantages when sudden changes in the mean flow occur, for surfaces 
with streamwise curvature, and in three-dimensional boundary layers. Both models show the Favre mass-weighted 
dependent variables account well for compressibility, even with rather large ap/3x, and also can account for 
modest effects of transverse curvature. 

Tnese models still require adjustments to increase their breadth of application. For example, the two- 
equation model needs an extra rate of strain added to the mixing-energy equation to account for the effect of 
streamwise curvature on a boundary layer. This ad hoc correction is very successful, practically, for 
boundary- layer calculations. In the Navier-Stokes form of the model, however, it has not yet been made to 
account for rapid turning within a flow. Nevertheless, the trail ing-edge example cited did not seem to need 
this correction. Both models are unable to completely relaminarlze an incompressible boundary layer in strong 
favorable pressure gradient. This is not a general failure of second-order models; another second-order 
closure model (Ref. 13) has been somewhat more successful in accounting for relaminarization than the models 
given here. Both models also require major changes in their surface boundary conditions to account for sur- 
face mass transfer or roughness. Finally, second-order models still require special treatment in regions 
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approaching irrotational flow. In conclusion, then, although the second-order closure models have shown a 
broader range of application than simple mixing-length methods, they are not universal and need further 
development to broaden their range of application, especially for boundary-layer flows that interact 
strongly with the surrounding irrotational flow. 

LARGE EDDY SIMULATION OF HOMOGENEOUS TURBULENCE 

The background and status of the techniques of large eddy simulation were recently reviewed by Ferziger 
and Leslie (Ref . l 44) . They gave particular attention to the methods for modelling the subgrid stresses 
expressed as <uiuj> in Eq. (14). To demonstrate the realism attained with the technique, the results of 
channel-flow computations were compared with data for mean-velocity profiles, Uj in Eq. (13), and mean 
moments, such as the pressure strain correlations. A more complete analysis of the simulation of channel 
flow appears in the Kim and Moin paper of this conference (Ref. 45). To supplement those paoers and to 
demonstrate the value of curbulence simulation to statistical turbulence modelling, this author will examine 
what can be learned regarding statistical Reynolds stress modelling from large eddy simulations of localized 
flow situations. 

The author is indebted to his colleague Dr. Robert Rogallo who generously provided the results of com- 
putations he is performing cn homogeneous turbulence that is experiencing decay, normal straining, or uniform 
shear straining. The Rogallo code has been described in Ref. 46, but has since been modified to accept uni- 
form shearing. The code has certain unique features. The turbulence is computed in a volume of fluid that is 
followed in time and is defined by coordinates that move with the assigned mean velocity. In this moving 
frame of reference, the turbulence is spatially homogeneous. The boundary conditions on the computational 
volume are treated as periodic in space, which permits use of full spectral methods in the computations. 

The code is efficient and accurate because particular care has been exercised to conserve energy and minimize 
aliasing. All variables are expressed in dimensionless form and related to physical quantities (subscript e) 
with scaling coefficients a and 6 as follows: 

Wave number or length 

Energy 

Kinematic viscosity 
Time 


It should be noted that the kinematic viscosity is treated as constant. 

In operating the Rogallo code, the turbulence is initially assigned an overall intensity with an arbi- 
trarily assigned three-dimensional spectral distribution. In addition, the mean strain rate and kinematic 
viscosity are also assigned. The turbulence is then oriented in phase space randomly while conserving mass. 
Because of the use of random phase, the components of turbulence velocity in each direction are uncorrelated 
so that no shear stress exists at time ■ 0. In the presence of a mean shear, the shear stresses develop in 
a short time and the computed results become independent of the particular random phase distribution that was 
used to start the calculations. The initially assigned spectrum also readjusts to be consistent with the 
assigned strain rate and kinematic viscosity, and the instantaneous turbulence intensity. 

The spectral range used in the calculations shown here has a ratio of the maximum to minimum wave number 
equal to 31. This ratio is established by (1) the storage capacity of the ILLIAC IV computer, which permits 
computations over volumes in phase space having 64 mesh points in each of three directions; and (2) the need 
for using two mesh spacings to define the minimum resolvable wavelength. Although this represents a very 
large number of computational mesh points, this spectral range is still inadequate to capture the range of 
wave numbers that is significant in a real turbulent flow, except for one at very small turbulence Reynolds 
numt?rs. Capturing the bulk of the significant eddy sizes and avoiding the use of a subgrid model is called 
an "honest" calculation. The Reynolds number appropriate to an "honest" calculation is an order of magnitude 
or more smaller than exists even in small scale laboratory experiments. If an "honest" calculation was to be 
compared with a low Reynolds number experiment, the physical output of the computations would be found from 
the calculations through Eqs. (50) after establishing a and B from the values of v and E used in the 
calculations and the v e and E e of the experiment. An alternative interpretation of these "honest" calcula- 
tions is to consider v used in the computations as an effective viscosity, which from Eqs. (13) and (14) is 
equivalent to the use of a constant eddy-viscosity subgrid model 

■ (v eff - v)(U 1,J + °1 ,J J (51) 

Emphasis is placed on the computation of the largest eddies in the flow, with the larger effective viscosity 
and a higher than real spectrum at the upper end of the wave numbers used to account for the dissipation that 
actually takes place at the wave numbers well beyond those in the computation. This approach presumes that 
the distorted spectrum at the upper end of computed wave numbers does not significantly alter the cascade of 
energy out of the energy-containing eddies at the low-wave-number end of the spectrum. In this approach, 
the scaling parameters o and 6 in Eqs. (5) can be established from comparing the calculated large eddy 
characteristics and corresponding quantities found in an experiment. 

To test the validity of this alternative interpretation of Rogallo's calculations, the computed results 
from several cases were compared with data obtained in the experiment by Harris et al. (Ref. 47) in which a 
rather complete set of turbulence measurements was made in a nearly homogeneous shear flow. When the com- 
puted macroscales and turbulence kinetic energy were matched to the corresponding experimental quantities in 
the region where the experiment reached an asymptotic behavior, the a and B needed to utilize Eqs. (50) 
were established. 

It was found that the best agreement between computation and experiment occurred when the effective 
viscosity in the computation was 15 times the molecular kinematic viscosity. A comparison of several computed 


6k , L e * B^L 
a“ *E 


t e = a 1 / 2 8*'t 


(50) 
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and measured statistical properties Is presented In Table 1. The computations show a kinetic energy dissipa- 
tion rate that Is about 20* higher than the rate Inferred from the measurements. This difference may not be 
significant however, because the energy dissipation rate Is difficult to measure accurately. It Is made up 
of direct measurements plus Inferences regarding local Isotropy of the small dissipative wave numbers. The 
excellent agreement of the mean velocity gradients Indicates the large turbulence structure is properly 
related to the mean flow In the computation. The remaining excellent comparisons between the Reynolds stress 
quantities Is strong evidence that the computation Is capturing the larger eddy structure that Is principally 
responsible for these quantities. This conclusion Is further supported by the comparison of the computed and 
measured two-point correlation coefficients for u. In the three Cartesian directions that are shown in 
Fig. 16. The general character of the experimental curves In each direction is represented very well by the 
calculations. Some of the weaknesses of the calculations are also demonstrated in this figure. The curve of 
R u (r,o,o) Is higher than that plotted from the data for the smaller separation distances. This is an indi- 
cation that the smaller eddies or high-wave-number eddies In the experiment are not well represented by the 
computations. This was expected In a computation with limited resolution, and one in which emphasis is on 
accurately computing the larger eddies. In addition, it Is noted that R n (r,o,o) has not vanished at 
rj/Lj = 4.9, which corresponds to 1/2 the length of each side of the computational volume. This suggests the 
computational volume used may have been too small and that the largest eddies could be sensitive to the 
periodic boundary conditions that were imposed. Even with these shortcomings, the remarkably good agreement 
between the experimental data and the computations encouraged the author to utilize the computations as a data 
base to examine some of the assumptions employed in statistical Reynolds stress modelling. 
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For a uniform homogeneous turbulent shear flow, the Reynolds stresses are 
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(52) 

(53) 

(54) 

(55) 


Closure of these equations requires expressing the correlation of pressure and rate of strain and the dissi- 
pation terms containing v in terms of the Reynolds stresses, themselves, and the mean flow. The pressure 
rate-of-strain model to be evaluated here Is represented by Eqs. (36) and (37). The dissipation terms in 
Eqs. (51) through (54) are usually replaced by the symbols tn, 622 , E 33 , and ei 2 . In terms of these quan- 
tities, the dissipation of the turbulence kinetic energy is given by 


e 


\ (< 


n 


+ £ ” + ‘ 33 * 


(56) 


It is usually assumed in modelling that the dissipation takes place at the smallest eddies and is, therefore, 
an isotropic phenomenon represented mathematically as 


£ 1j*f £S 1j (57) 

To learn if the computations are consistent with the assumption of isotropic dissipation in a shear 
* flow, the computed values of E-jj corresponding to each Reynolds stress are plotted against a measure of the 
anisotropy, u^uj/e - 2/3 5 1 j , in J Fig 17. Ihe alues corresponding to the different Reynolds stresses for a 
range of turbulence Reynolds numbers gener-ill;, lie along a straight line with a slope of about 0.7. If the 
dissipation had been isotropic, in the coordinates of the figure, these points would have been located on the 
axis or at an ordinate equal to zero. The computations show the dissipation to be anisotropic and require 
that Eq. (57) be modified to 


E 1j 



+ 0.7 



(58) 


for the shear flow. What is normally termed dissipation in a large eddy simulation actually represents the 
component energies that are cascaded toward the high wave number end of the calculation to be then drained 
from the calculation by the subgrid model. It is no surprise then that the cascade process reflects the 
anisotropy of the larger eddies. The emphasis on t.he behavior of the large eddies possibly is an advantage: 
it may be just what is required in Reynolds stress modelling, which also addresses the behavior of the 
larger eddies. 

Figure 18 shows the components of dissipation as functions of anisotropy for the case of a turbulent 
flow relaxing after it had been instantaneously distorted with normal strains in the x 2 and x 3 directions. 
In this case, the computed results generally lie along the coordinate axis and the dissipation is approxi- 
mately isotropic even though the flow Itself is sti‘1 anisotropic. 

The homogeneous flow relaxing after Instantaneous distortion by normal strains is also an excellent 
case for examining the first term in the pressure rate of strain relationship represented by Eq. (36). The 
terms in Eq. (36) preceded by a, 6 , and y are identically zero in the absence of continuing mean strain. 

The C t were evaluated from the computed turbulence moments for a case distorted in all three directions, 
it was found that Cj was reasonably insensitive to the direction of the component considered. The Cj 
(based on U 2 ) Is plotted as a function of the turbulence Reynolds number e 2 /Ev e ff in Fig. 19. The numbers 
adjacent to the symbols represent the magnitude of the largest anisotropy in the three components. Three 
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different values of v e ff were used In the calculations to extend the range of the turbulent Reynolds number. 
It Is observed that If the largest component of anisotropy Is less than 0.25, that Is, If 



< 0.25 


(59) 


the values of Cj collapse onto a single curve. For those values of Isotropy where the points depart from 
a single curve, the linear form o, the Potta term expressed in anisotropy can be considered to have failed. 
It Is also noted that for values of anisotropy satisfying (59), the so-called "constant" Cj still has a 
strong Reynolds number dependence. 


The value of Cj was also evaluated from the homogeneous shear flow computation of the pressure rate- 
of-strain correlation and Eqs. (36) and (37). To do this most simply, the interrelationships between a, 

B, and y derived by Launder et al . (Ref. 14) following Rotta's suggestion (Ref. 8) were utilized. The Cj 
formed from these computations is plotted as a function of turbulence Reynolds number in Fig. 20. Again, 
the maximum anisotropy is indicated at the plotted points. For comparison, the C ; from the normally 
strained flow, with the value of maximum anisotropy characteristic of shear flow, are also plotted on the 
figure. Shearing the flow tends to increase the turbulence Reynolds number so that the regions of the two 
flows do not overlap. The same magnitude of anisotropy is used for both flows to account for a similar 
departure from the linear Rotta form in each. It is observed that the line segments associated with the 
different types of flow fields do not appear to form a common curve. This would imply that the pressure 
rate-of-strain model represented by Eqs. (36) and (37) is not as universal as is suggested by its tensor 
form. It is most Interesting, however, to note that if the Cj is considered to be the coefficient of the 
sum of the pressure rate-of-strain correlation and the anisotropic dissipation, the dashed curves on Fig. 20, 
the two flow fields tend to produce a common curve. It appears that only the sum of the dissipation and 
pressure rate-of-strain can be modelled universally. This observation is consistent with the theoretical 
approach adopted by Lumley and Newman (Ref. 48). It should be noted, again, that the anisotropy on the 
figure is outside the region of applicability of the Rotta relationship and that shear flows with lesser 
anisotropy would require larger values of Cj for the pressure rate of strain contribution. 

The values of Cj, a, B, and y that conform to the high end of the turbulence Reynolds number in these 
calculations are shown in Table 2 along with those used in the Launder et al . (Ref. 14) model and the RSE 
model described here. The agreement between the computations of Launder et al . and those of Rogallo is 
excellent and is primarily due to the ability of the computations to yield good values of the Reynolds 
stresses in equilibrium shear flow (Table 1). The requirement of i M in the RSE model apparently requires 
considerable compensation in all the other terms to result in the proper Reynolds stress ratios for uniform 
shear flows. 


Concluding Remarks Regarding Large Eddy Simulations 

Although this demonstration of the use of large-eddy simulations for guiding Reynolds stress modelling 
has been limited to homogeneous flows, its utility as a research tool shows great promise. Although the 
procedure is quite costly in terms of computer time (the calculations shown here require about 1.5 hr 
ILLIAC IV time per test case), the potential gains in computer technology (Ref. 1) should make large eddy 
simulations a research too: that will be available to most research laboratories in a decade or so. It is 
this author's belief that research activities Involving coordinated theory, experimentation, and computer 
simulations will in the reasonably near future not only bring about a much clearer understanding of the 
physics of turbulence, but may even permit the development of predictive engineering methods for flow fields 
of technological interest. 
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TABLE 1. COMPARISON OF ROGALLO'S SHEAR FLOW CALCULATIONS 
WITH THE EXPERIMENTAL DATA OF HARRIS, GRAHAM, AND C0RRSIN“ 


Quantity 

Experiment 

Computation 

Dissipation rate, cm 2 sec" 3 
Mean- velocity gradient, Uj 2 sec -1 

3.28E+04 

3.92E+04 

44.0 

45.3 

Angle of principal stressed, deg 

-22.3 

-22.6 

Ratio of principal stresses 

4.06 

5.24 

UjUj/e 

1.00 

1.01 

u 2 u 2 /e 

0.40 

0.36 

UjUj/e 

0.60 

0.63 

-uiu 2 /e 

0.30 

0.33 


“Scaling established by matching turbulence kinetic energy 
and streamwise macroscales. Turbulence model v e ff/v ■ 15. 


TABLE 2. PRESSURE RATE OF STRAIN CORRELATION 
MODELLING COEFFICIENTS 



Launder, Reece, and 

Wilcox and 

Rogallo' s 


Rodi model 

Rubesin model 

computations 

a 

1.5“ 

4.5“ 

1.5“ 

0.76 

0.5 

0.78 

B 

0.11 

0.5 

0.23 

Y 

0.36 

1.33 

0.55 


“includes effects of anisotropic dissipation. 
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Fig. 1. Comparison of computed and measured dis- 
tortion parameter for the Tucker- Reynolds plane 
strain flow. 


c t 




(a) Skin- friction distribution. 

(b) Shape-factor distribution. 



u T y/r 


Fig. 3. Comparison of computed and measured velocity 
profiles for the Bradshaw adverse-pressure-gradient 
flow; x * 2.1 m. 


■ BRADSHAW DATA 

RSE MODEL 

TWO-EQUATION MODEL 



Fig. 2. Comparison of computed and measured skin 
friction and shape factor for the Bradshaw adverse- 
pressure-gradient flow. 


Fig. 4. Mean field profiles for the Bradshaw adverse- 
pressure gradient flow; x * 2.1 m. 
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Fig. 5. Effect of transverse surface curvature on 
skin friction. 




(a) Skin- friction distribution. 

(b) Shape-factor distribution. 

Fig. 6. Comparison of computed and measured skin 
friction and shape factor for flow over a convex 
wall with adverse pressure gradient. 


Fig. 7. Effect of surface cooling on the skin fric- 
tion on a flat plate. 


O DATA OF LEWIS, GRAN, & KUBOTA 
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Fig. 8. Skin-friction coefficient distribution in a 
region of large pressure gradient at M«, * 4. 










11-21 


NOZZLE 


» TEST SECTION 


! 24. 77 


X DIFFUSER 
cu/v*k 

4 GENERATOR 

4 


o EXPERIMENT 
TWO EQUATION MODEL 


RANKINE-HUUONIOT . 



Fig. 11. Comparison of computations and surface 
measurements of pressure and skin friction In 
region of shock wave, boundary- layer Interaction at 
low supersonic speeds. 
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Fig. 12. Comparison of computations and surface 
measurements of pressure and skin friction in the 
region of a compression comer at supersonic 
speeds. 


EXPERIMENT 
o o ■ 20 

o o-16 

ZERO-EQUATION MODEL 

ONE-EQUATION MODEL 

TWO-EQUATION JONES-LAUNDER MODEL 

TWG-EOUATION WILCOX-RUBESIN MODEL 

1.6 



Fig. 13. Comparison of computations and measurements 
of the effect of Reynolds number on the extent of 
upstream pressure Influence of a supersonic compres- 
sion corner, M * 2.8. 

M„ - 0.7, Re L - 36.6 X 10®, o = 0 deg 

O EXPERIMENT 

ALGEBRAIC MODEL 
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Fig. 14. Comparison of computed and measured mean- 
velocity profiles upstream and downstream of sharp 
trailing edge with a 12.5° wedge angle. 
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Fig. 15. Ratio of the transverse and longitudinal eddy viscosities In the region of a sudden application of 

transverse shear. 




Fig. 16. Two-point correlation coefficients of Fig. 17. Components of dissipation within a homoqe 

u ] along the three Cartesian axes. neous shear flow. 
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Fig. 18. Components ov' didsipatlon In a relaxing 
homogeneous turbulence after sudden normal 
distortion. 


Fig. 19. Behavior of the coefficient in the Rotta 
pressure rate of strain correlation term in a relax- 
ing homogeneous flow after sudden distortion by 
normal forces. Numbers near symbols represent the 
value of the largest component of anisotropy present. 



Fig. 20. Reconciliation of modeling the Rotta term in normally and shear strained flows. Numbers near the 
symbols represent the value of the largest component of anisotropy present. 


